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Abstract 

One of the aspects currently holding back commercial scale deploy¬ 
ment of carbon capture and storage (CCS) is an accurate understanding 
of the thermodynamic behaviour of carbon dioxide and relevant impurities 
during the pipeline transport stage. In this article we develop a general 
framework for deriving pressure-explicit EoS for impure CO 2 . This flexi¬ 
ble framework facilitates ongoing development of custom EoS in response 
to new data and computational applications. We use our method to gen¬ 
eralise a recent EoS for pure CO 2 [Demetriades et al. Proc IMechE Part 
E, 227 (2013) pp. 117] to binary mixtures with N 2 , O 2 and H 2 , obtaining 
model parameters by fitting to experiments made under conditions rele¬ 
vant to CCS-pipeline transport. Our model pertains to pressures up to 
16MPa and temperatures between 273K and the critical temperature of 
pure CO 2 . In this region, we achieve close agreement with experimental 
data. When compared to the GERG EoS, our EoS has a comparable level 
of agreement with CO 2 -N 2 VLE experiments and demonstrably superior 
agreement with the O 2 and H 2 VLE data. Finally, we discuss future op¬ 
tions to improve the calibration of EoS and to deal with the sparsity of 
data for some impurities. 
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1 Introduction and Background 

1.1 Carbon capture and storage 

Carbon capture and storage is a crucial technology in the international efforts 
to meet carbon dioxide emission targets [1]. Capturing CO2 from industrial 
sources can lead to a 90% reduction in emissions. However, no gas separa¬ 
tion process is 100% efficient, and as a result the CO2 generated from power 
generation or by industry can contain a number of different impurities, depend¬ 
ing on its source. These impurities can, depending on their composition and 
concentration, greatly influence the physical properties of the fluid compared 
to pure CO2. Impurities have important design, safety and cost implications 
for the compression and transport of CO2 and its storage location, for example 
geological sequestration. This research is designed to tackle one of the key tech¬ 
nical challenges facing the development of commercially viable CO2 transport 
networks: modelling physical behaviour of impure CO2, under the conditions 
typically found in carbon capture from power stations, and in high-pressure 
(liquid phase) and low-pressure (gas phase) pipelines. Accurate modelling of 
the physical properties of CO2 mixtures is essential for the design and opera¬ 
tion of compression and transport systems for CO2. In particular the variation 
of fluid density and phase-behaviour with temperature and pressure is key to 
many CCS processes. For example, pipeline transport of CO2 is only viable if 
the fluid remains in the homogenous phase. 

There has been recent work to define the expected operating conditions 
for CCS pipelines [HIS]. The most efficient way of transporting CO2 is in the 
homogenous phase, at pressures in the vicinity of its critical point. For the 
transport temperature, the upper temperature will be set by the compressor 
discharge temperature and the temperature limits of the pipeline coating and 
the lower temperature will correspond to the winter ground temperature of the 
surrounding soil[3]. Expected impurity levels are about < 4%, with N2 , O2 and 
H2 being key impurities [llllili. This range of pressure, temperature and 
impurity level define pipeline operating conditions and provide a target window 
for CCS-oriented modelling. However, CCS-relevant models should aspire to 
model a wider range of conditions, particularly for the impurity level, for the 
following reasons. Coexistence leads to the formation of a vapour phase that 
is considerably richer in impurity than the overal mixture. Upset conditions, 
in which a greater concentration of impurity is accidentally introduced into the 
CO2 stream, must be understood and mitigated for. Finally, an effective way to 
ensure physical robustness of the model is to test for a wider range of impurity 
conditions. 

1.2 CO 2 modelling in CCS 

In models of the CCS process, the fluid behaviour of CO2 mixtures is typically 
predicted by an equation of state (EoS). EoS vary in their mathematical form, 
accuracy, region of validity and computational complexity. Because different 
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applications have different requirements there is no single EoS that is ideal 
for all applications. In particular, there is a balance between mathematical 
simplicity and accuracy of prediction. To optimise their accuracy, EoS need 
to be calibrated by fitting their parameters to experimental measurements on 
CO2 mixtures. However, new measurements become available very frequently, 
offering the opportunity to improve the models. Thus, there is an ongoing need 
to regularly rederive, refine and reparameterise EoS. Currently, the inability to 
rapidly assimilate ongoing measurements into suitable EoS delays or prevents 
knowledge gained from experiments from being applied in CCS modelling. 

A pressure-explicit EoS is an expression for a fluid’s pressure P*, as a func¬ 
tion of molar volume v* and temperature T*. Usually, the terms in the model 
that describe the deviation from ideal gas behaviour are empirically postulated. 
A widely used example is the Peng Robinson equation [ 7 ]. Such EoS gives 
an explicit prediction of the pressure-volume curve for homogeneous fluids and 
can also predict the coexistence behaviour. For thermodynamic co-existence, 
the coexisting phases must have matched fugacity for each chemical species, 
which is equivalent to matching the chemical potential. This translates into a 
constraint involving the integral of P* over volume, at constant temperature. 
Thus, the coexistence behaviour can be predicted by numerically searching for 
two volumes that obey both the EoS and thermodynamic coexistence require¬ 
ment. Typically, EoS contain empirical parameters that are estimated by fitting 
the EoS to experimental data for the pure material. These parameter are gen¬ 
eralised to mixtures through a set of empirical mixing rules that extend the EoS 
to multi-component mixtures. These mixing rules also allow calculation of the 
phase behaviour, through a generalised expression for the fugacity. The pure 
component parameters and the mixing rules may vary with temperature. 

There is an alternative method to formulate EoS, in which the volume and 
temperature-dependence of the Helmholtz free energy, rather than the pressure, 
is postulated^ | 9 ]. Similarly to pressure-explicit EoS, these models contain 
empirical terms that described the deviation from ideal gas behaviour. The 
two formulations are, however, mathematically equivalent as integration of a 
pressure-explicit EoS leads to an expression for the Helmholtz free energy. In 
this work, we use the pressure-explicit formulation, as measurements of pres¬ 
sure can be directly visualised and so postulating empirical terms for P*{v*) is 
physically more intuitive than for the free energy. 

There is considerable uncertainly over which EoS is most appropriate for 
CCS modelling. Although the Peng-Robinson model is mathematically sim¬ 
ple and numerically cheap, it was derived to compute separation of mixtures 
for the natural gas industry, where CO2 is a minor additive. Thus, its pa¬ 
rameters are optimised to the coexistence behaviour, particularly for the gas 
phase. Agreement with density measurements for pure CO2 around the criti¬ 
cal pressure, which is key to CCS modelling, is unacceptably poor. Moreover, 
the mixing rules are not optimised to CCS-relevant mixtures. Variants of the 
Peng-Robinson model usually suffer the same limitations. 

For pure CO2, the Span-Wagner EoS [ 8 ] covers from the triple-point tem¬ 
perature up to very high pressures and temperatures with very high accuracy. 
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Furthermore an EoS by YokozekilTU] captures solid-liquid coexistence of pure 
CO2. Also for pure CO2, there is a recent composite EoS[TT] that combines 
the Peng-Robinson model for the gas phase and accurate tabulated data for 
the solid and liquid phases [HI HU US]- There have also been studies comparing 
different pure CO2 EoS for CCS applications mm- There are complex and ac¬ 
curate EoS for CO2 mixtures, including the SAET[TH], PC-SAET[TH], GERCjS] 
and EOS-CG[T 71 [TH] models. The SAFT, PC-SAFT and EOS-CG models have 
recently been compared to some CCS-relevant measurements [lH I 18 j . How¬ 
ever, they have not been compared to recently emerged VLE data for CO2 -H2 
mixtures |201 [ 5 ], Furthermore, the mathematical complexities of these models 
preclude their use in some CCS applications. 

There is clearly considerable scope to derive new EoS, specifically for CO2 
mixtures, retaining the simplicity of Peng-Robinson type equations but with 
improved quantitative performance for CCS applications. Indeed, a simplified 
version of the Span-Wagner model for pure CO2 has been produced jH]. How¬ 
ever, the measurements against which a CCS-focussed EoS should be calibrated 
continually evolve as CCS-oriented projects deliver new data to complement the 
literature data. Additionally, every user’s requirements differ, depending on the 
accuracy demanded, the region of interest, and the computational complexity 
that can be tolerated. Thus, there is no single “silver-bullet” EoS, suitable for 
all CCS applications and all users. Instead, each application’s requirements can 
only be met by a bespoke EoS. There is a need for a flexible, general framework 
to derive and parameterise new EoS in response to emerging measurements and 
computational applications. 

1.3 Aims of this work 

In this work we introduce a general and flexible framework for deriving new 
pressure-explicit EoS. We derive an expression for the mixture fugacity for an 
arbitrary EoS with arbitrary mixing rules. This allows an EoS to be generalised 
to mixtures, for any choice of mixing rules, without needing to recompute the 
fugacity integral. Therefore, our approach allows convenient modification of the 
form of EoS and mixing rules, in response to new measurement or new physical 
insights. We demonstrate the flexibility and effectiveness our approach by gen¬ 
eralising a recent pressure-explicit EoS for pure CO2 to mixtures with N2, O2 
and H2 and then calibrate this model against data for CCS-relevant mixtures, 
focussing on pressure-volume behaviour and coexistence. Our approach is read¬ 
ily generalisable to other common CCS impurities, such as argon and methane 
and can also fit to other thermodynamic variables such as enthalpy and heat 
capacity. 
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2 A general expression for the mixture fugacity 
in pressure-explicit equations of state 


In this section we consider a general pressure-explicit EoS and derive results 
for arbitrary EoS and mixing rules. To compute coexistence we require an ex¬ 
pression for the fluid fugacity. For pure fluids this can be obtained by directly 
integrating the EoS, using eqn (15). Whether this integral can be performed 
explicitly depends on the expression chosen for the EoS. Furthermore, for the 
mixture fugacity, one must compute a more complicated integral, involving the 
mixing rules (eqn (16l). Computing this mixture integral directly is not mathe¬ 
matically convenient for all but the simplest EoS. In this section we show that, 
provided the pure fluid integral can be performed, then a closed formed expres¬ 
sion for the mixture fugacity integral can always be found. Furthermore, we 
derive a general expression for the mixture fugacity in terms of derivatives of 
the mixing rule. 

We begin with an EoS in the form P* = P*{v*,0*,T*), where P* is the 
pressure in Pa, v* is the molar volume in mol/m^, T* is temperature in K and 
9* is a vector of model parameters. Here, stars denote dimensional quantities. 
The fugacity co-efficient for the pure fluid, Ini^, is given by eqn (15) in|^ The 
EoS is generalised to mixtures by allowing the parameter 0 to depend on the 
fluid composition, x, to be specified via a mixing rule. Here Xi denotes the mol 
fraction of species i of the mixture. Leaving the functional forms of P* and 0* 
unspecihed, we show in that the fugacity coefficient of species z in a mixture 
is given by 


In (j)i{v*,x) 


Np 

ln(()-b^ 

i=i 


dP 



Nsp 



(1) 


where F = Int/) — {{Z — 1) — InZ). We see that, if the fugacity integral can 
be performed for the pure fluid, then eqn ([^ ensures that the fugacity inte¬ 
gral for mixtures can also be written as a closed form expression. Finally, we 
note that eqn 0 is invariant under a constant scaling of any of the model 
parameters. This means that the expression is completely unchanged by any 
non-dimensionalisation that merely scales the model parameters by a constant. 


3 A new Equation of State for impure CO 2 


3.1 Pure CO 2 

We begin with a recently proposed EoS for pure CO 2 by Demetriades et al. [52] . 
The dimensional EoS is 

RT* 6*^ d*^ 

P(v* T*-0*) = — _°_ ^ _h 

’ ’ V* + a* v*'^ c*^ v*^ + e*^ 
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where R is the ideal gas contant and 0* = {a*,...,g*) is the vector of model 
parameters. We perform the following non-dimensionalisation, for the physical 
variables 


P* = P*P, T* 
and for the model parameters 


r*r, 



(3) 
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Here P* and T* are the critical pressure and temperature of pure CO 2 , respec¬ 
tively. After non-dimensionalisation our EoS becomes 


P(n,T) = 


T 


V- g 


V + a 

and the fugacity integral for pure CO 2 (eqn ( [I^ ) becomes 

which can be evaluated by substituting in eqn ([^ 
ln0(ti) = In 
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f 

hT{g - vf 

(7) 


We repeated the fitting procedure for pure CO 2 from Demetriades et a/.[^ 
to obtain slightly improved agreement with the pure CO 2 data over our previous 
article. The resulting variation of model parameters with temperature is given 

by 


aco^T) 
bco2 (T) 
CC02{P) 
dco2{P) 

6002 ( 7 ") 

fco2{T) 

5002 iP) 


\T - 1|° ®2®207('33 9202^12^ _ _ 8.10461|T - 1| -h 0.805812) -I- 0.2712941, 

\T - 1|° '‘°®254^_j^ 3 52Qg|2. _ ^|2 4_48534|r _ i| _ 0.295229) -f 0.3326169, 

\T - 1|°-®1®^®®(-3.77054|T - 1\^ + 1.72673|r - 1| - 0.478733) -f 0.238762, 

\T - l|^-2^°®®(0.000634507|r - if - 8.327888244017052E-6|T - if 
-0.0000382867|T - 1| -h 4.661593764290955E-6) - 0.000374407355, 

0.780746514, (8) 

\T - 1|°-^^2269('q 2io 429|T - Ip - 0.199813|r -1\+ 0.0528131) -f 0.0787701, 

\T - l|°-^^®^^^(-0.185594|T - Ip 0.0931741|T - 1| - 0.0510056) -f 0.074028115. 
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Thus each parameter was fully defined and this completed the model in the case 
of pure CO 2 . We note that the parameter variations given here differ slightly 
from those proposed in our published work |22]. The model proposed here 
has a better accuracy than that of our prior version. We now fix these model 
parameters for the remainder of this article. 

3.2 Generalisation to impure CO 2 

We generalise our model to impure CO 2 by allowing the model parameters to 
depend on the impurity concentration, a^imp- We opted to impose a linear mixing 
rule for each of the model parameters in equation ([^. For parameter a, the 
mixing rule for a binary mixture of CO 2 and a single impurity (imp) is: 

a(T) — a:c;02^C02 (^) T ^impaimp(T), (9) 

where aco 2 {"r)is defined in equation (|^ for pure CO 2 and ainip(T) is the analo¬ 
gous quantity for the impurity, which we need to specify. Replacing a in eqn 
with any of the other model parameters provides the mixing rule for that pa¬ 
rameter. We note at this stage that for a binary system 

^imp 1 ^C02 ■ 

From here on we use x to denote the impurity fraction. 

3.3 Fitting strategy for impurity data 

Our htting strategy for the impurity model is similar to our methodology for 
pure CO 2 [22] . in that we fit model parameters at different temperatures by using 
simulated annealing, which searches for global minima and has the ability to 
escape local minima. The correct coexistence behaviour is ensured by imposing 
a term in the minimisation that penalises differences in the predicted fugacity at 
the experimentally determined coexistence points (see eqn This approach 

requires vapour-liquid equilibrium (VLB) measurements, for both the coexisting 
mol fractions and molar volumes at each temperature to be fitted. 

4 Carbon Dioxide—Nitrogen Binary Mixture 

4.1 Data Availability 

We begin this section by noting the literature data available to fit the CO 2 -N 2 
mixture, which is summarised in TableAs specihed in above, we require VLB 
data, including coexisting volumes, ideally with homogeneous phase density data 
all at the same temperature. There are only two relevant temperatures in the 
literature where this full spread of data exist and, in some cases, the fraction 
of N 2 used for the homogenous phase measurements is outside the range of 
relevance to CCS. In order to find the pure N 2 parameter values, we fitted our 
model to the CO 2 -N 2 mixture data at 273.15K and 288.15K, choosing these 
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Temperature (K) 

Homogeneous Density 

VLE Data 

260 

m, m, m 


265 

m, m 


270 

m, m 

Ea^ I3a^ 

273.15 

1^. na, M 

EH, [^1, [Ml|36],E7]i , 



EH]i, Ea^ Ila^ 

275 

m, m, m, na, m 


280 

m, m 


285 

m, m. EH, m, Ea 


288.15 

m 

m, I3a^ mi^ 

288.706 

Ea 


290 

m, m 


293.3 

na 

nut 

295 

m, m 


298.15 

sa 

El]i 

300 

m, m. m, m. m 



m, Ea, na, m 


301.3 


m 

303.3 

m 

m 

^ means the VLE data in this source do not contain volumes 


Table 1: A summary of available data for CO 2 -N 2 mixtures. Experimental 
data used in our parameter fitting are in boldface. 


temperature as they both had measurements of the coexisting volumes and 
homogenous phase. A summary of the data used for fitting is as follows. For 
273.15K we used VLE data by Muirbrook et and homogenous phase 

data by Aral et ol.[32] for xm^ in the range 7.1 — 25.4%. For 288.15K we used 
VLE data by Aral et aL[32] and homogenous phase data by the same authors 
for in the range 6.0 — 10.3%. 

4.2 Fitting the CO 2 —N 2 Binary Model Parameters 

We proceeded first by determining the pure N 2 parameters onj to at the 
available temperature points by a similar series of optimisations as used in the 
pure CO 2 case. Specifically, we minimised an error function, which in the case 
of a binary mixture, was comprised of two contributions: from the homogeneous 
measurements 


A'Hom — ^ ^ 


T 1 — P' 

( VdATA’^ ’-^dataI m 


(0 

DATA 


i=l 




Pt 


1 ) 


DATA 


( 11 ) 


where is the number of homogeneous measurements at this temperature, 
Upj^TA ^od a:DAxA denote the volume and N 2 composition, respectively, of the 











0 

ao 

Cti 

aN2 

-2.94804 

3.14877 

0N2 

0.616547 

-0.654773 

CN 2 

-4.83602 

5.83048 

C?N2 

4.69708 

-5.23021 

eN2 

6.23314 

-6.56762 

/n2 

-11.6889 

13.4184 

ffN2 

15.3391 

-17.2528 


Table 2: Coefficients describing the temperature dependence of the N 2 model 
parameters in equation (14). 


ith homogenous data point and P{v,T,x) is the model’s predicted pressure. 
The second contribution comes from the VLB measurements 


^^VLE = 


(P(4i)p, T, T, xg) - 


E 

i=i 

+ (ln^C02(4iL4il) - ln^C02K?n\4q) + ln((l - 4ip)/(l “ 4u))) 


+ (ln(/.N2(w™p,4aU - ln0N2 


( 12 ) 


where is the number of VLB measurements at this temperature, 

and a^l^p/iiq denote the coexisting liquid/vapour volume and N 2 composition, 

respectively, of the jth VLB measurement, Pvip is the jth vapour pressure 
measurement and In ^ is the model’s predicted fugacity. The two contributions 
were combined into a total error function 


ETot.l = WEnom + (1 - W^)^ivLE, (13) 

where VT is a parameter that determines the weighting of the fitting between 
homogenous and VLB data. We chose W = 0.1, to balance fitting between the 
VLB and homogenous phase data. This generated values for each of the seven 
pure N 2 parameters at both of the temperature points. 


4.3 Variation of Model Parameters with Temperature 

We had determined values for each parameter at each of the two temperatures. 
To obtain an expression for the temperature dependence of each parameter we 
converted our parameter values into a linear dependence on temperature, 

9n^[T) = ao +aiT, (14) 

where 0 is a model parameter. The values of ao and ai we obtained for each of 
the seven N 2 parameters are in table 
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4.4 Performance of the Model for Binary Mixtures of CO 2 and 
N2 



Figure 1: Comparison of experiments and our model for the VLE coexisting mole 
fraction for CO 2 -N 2 mixtures. Experimental data at 273.15K from Muirbrook 
et at[35l[36] and at 288.15K from Aral et aZ.[32]. Also shown are the results 
from the GERG EoS [9]. 


With the model thus fully defined for the binary system involving CO 2 and 
N 2 , we were able to assess its performance. We did this by comparing the com¬ 
pleted model to the data was used to generate the fitted parameters. We com¬ 
puted the model’s coexistence predictions, at a given temperature, as follows. 
We began at the pure CO 2 coexistence pressure, where the coexisting volumes 
are known and the coexisting mol fractions are both 100% CO 2 . We then in¬ 
creased the pressure by a small increment and searched for the mol fraction 
and volume of the liquid and vapour that satisfy the following four coexistence 
conditions: the liquid phase pressure equals the prescribed pressure; the vapour 
phase pressure matches the prescribed pressure; the CO 2 fugacity of the liquid 
matches that of the vapour; and the N 2 fugacity of the liquid matches that of 
the vapour. We used a numerical non-linear root-finding algorithm to locate the 
co-existence point. Once the coexistence root was found we increased slightly 
the pressure and repeated the root finding, using the previous root as an initial 
guess for the new root. We repeated this process, moving up in small pressure 
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Coexisting volumes 
Homogeneous phase (Pure CO 2 ) 
Homogeneous phase (7.1% N^) 
Homogeneous phase (12.8% N^) 
Homogeneous phase (25.4% N 2 ) 



0.0001 


0.0002 

V* / 


0.0003 

(m^/mol) 


0.0004 


0.0005 


Figure 2: Comparison of experiments and our model for the pressure-volume 
behaviour of CO 2 -N 2 mixtures at 273.15K in both coexistence and the homoge¬ 
nous phase. Experimental data from Muirbrook et aL |35l 155] (coexistence) and 
from Arai et aL[32] (homogenous phase). 


increments, until the critical point was reached. 

A comparison between the experiments and our model for the coexisting mol 
fraction is shown in figure Also shown are the results from the GERG EoS. 
At both temperatures for both models the agreement is very good. Our model 
performs similarly to the GERG, except at pressures approaching the critical 
point for 288.15K. Here, our model fails slightly to capture accurately the ap¬ 
proach to the critical point, which is most clearly seen by the rapid convergence 
of the liquid and vapour mol fraction as the experiments approach the critical 
pressure, a feature that the GERG captures. For coexisting volumes ( figures 
and our model performs very well, except for close to the critical point at 
288.15K, where our model slightly over-predicts the pressure. Our model cap¬ 
tures the mixtures’ pressure-volume data for the homogenous phase (figures [^, 
except at high pressures at 288.15K where the model’s results are slightly closer 
to the pure CO 2 behaviour than the experiments. Finally, we also highlight that 
the best performance of our model occurred for the liquid coexistence proper¬ 
ties. This will ensure that the model accurately predicts the liquid saturation 
line, which is of particular relevance to CCS pipelines as it is the edge of the 


11 











V* / (m^/mol) 


Figure 3: Comparison of experiments and model predictions for the pressure- 
volume behaviour of CO 2 -N 2 mixtures at 288.15K in both coexistence and the 
homogenous phase. Experimental data from Aral et al.\62\. 


homogeneous phase region and, ultimately, it defines the minimum safe pipeline 
operating pressure for a given temperature and overall impurity composition. 

Overall, this mixture shows the promise of our methodology. In the areas 
where our current model fails a little to fully capture the data, namely close to 
the critical pressure, the flexibility and generality of our approach will enable 
future work to readily introduce refined mixing rules or terms in the original 
EoS to obtain further improved agreement. Conversely, our framework will also 
enable the complexity of the EoS to be reduced for applications that demand a 
simple EoS but can compromise on the closeness to experiments. 


5 Carbon Dioxide—Oxygen Binary Mixture 

5.1 Data Availability 

The availability of data for CO 2 -O 2 is more limited than for CO 2 -N 2 . A 
tabulated summary of the relevant available data is given in Table There is a 
single temperature at which the full array of VLE information was presented, at 
273.15K [^ES]. The highest temperature at which any VLE data is presented 
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Temp. 

Homogeneous 

Density 

VLE Data 

263.15K 

273.15K 

283.15K 
302.22K 

m 

[H] t, [351136], mi t 
nut 

^ these VLE data do not contain volumes 


Table 3: A summary of available data for CO 2 -O 2 binary mixtures. Experi¬ 
mental data used in our parameter fitting are in boldface. 


0 

OiQ 

Oi\ 

aoa 

-9.3985 

10.4231 

bo2 

4.62523 

-5.15272 

^02 

-1.50723 

1.9847 

do2 

12.838 

-14.292 

eo2 

5.99309 

-6.18695 

/02 

9.29129 

-9.9675 

0 

to 

-7.61445 

8.30386 


Table 4: Coefficients describing the temperature dependence of the O 2 model 
parameters in eqn (14). 


is 283.15K, and there is no temperature at which both VLE and density data 
were given. To enable fitting to a second temperature we developed a volume 
estimation technique for the co-existing volumes, as detailed in to fill-in 
where experimental measurements were unavailable. A summary of the data 
used for fitting is as follows. For 273.15K we used VLE data by Muirbrook et 
al. [3S1I3S]. For 283.15K we used VLE data from Fredenslund et a7[5T], with 
estimated coexisting volumes. At both temperatures we also used pressure- 
volume data for pure O 2 from the NIST website [K2] . 


5.2 Fitting the CO 2 —O 2 Binary Model Parameters 

To fit the CO 2 -O 2 parameters, we proceeded in exactly the same way as for 
fitting the impurity parameters in the previous section. We used the error 
function in eqn (13) and fitted the seven impurity parameters via simulated 
annealing, treating each temperature separately. We used W = 0.1, to balance 
fitting between the VLE and homogenous phase data. To obtain an expression 
for the temperature dependence of each parameter we converted our parameter 
values into a linear dependence on temperature (eqn (14|). The coefficients of 
this linear expression for each of the parameters are listed in table 
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Figure 4: Comparison of experiments and model predictions for the VLE 
coexisting mole fraction for CO 2 -O 2 mixtures. Experimental data at 273.15K 
from Muirbrook et a7[3Sl[35]; and at 283.15K from Fredenslund et aZ.[^. Also 
shown are the results from the GERG EoS [5] . 


5.3 Performance of the Full Model for Binary Mixtures of 
CO 2 and O 2 

Similarly, to GO 2 -N 2 mixtures we computed the VLE correlations of the model 
and compared to the data used for fitting. A comparison between the exper¬ 
iments and our model for the coexisting mol fraction at both temperatures is 
shown in figure]^ Also shown are the results from the GERG EoS. At 273.15K 
our model captures well the data, with a very slight discrepancy very close to 
the critical point. In contrast, the GERG dramatically over-predicts the critical 
point pressure and has noticeable disagreement with the liquid mol fraction that 
persists over most of the pressure range. There is a similar picture for the GERG 
at 283.15K, whereas our model under-predicts somewhat the critical point at 
this temperature. Overall our model produces consistently better agreement 
than the GERG. We note the SAFT and PC-SAFT EoS have been compared 
to these data at 288.15K[Tn] and both over-predict the critical point pressure to 
a simular extent to the GERG. Finally the EOS-GG model has also been com¬ 
pared with these experiments[T5]. The EoS-CG over-predicts the critical point 
pressure by ~ 20% (much less than the GERG and SAFT models) and has bet- 
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Figure 5: Comparison of experiments and our model for the coexisting volumes 
for CO 2 -O 2 mixtures at 273.15K. Experimental data from Muirbrook et al . [351 

m- 


ter agreement to the liquid mol fraction experiments than our model. Overall 
our model has superior agreement than the GERG and the SAFT models and 
comprable agreement to the EOS-CG model. 

Eigurej^ compares our model with measurements of the coexisting volumes 
for CO 2 -O 2 mixtures at 273.15K. Agreement is very close throughout. 


6 Carbon Dioxide—Hydrogen Binary Mixture 

Coexisting mol fractions for CO 2 -H 2 mixtures under conditions relevant to CCS 
have recently been measured by Eandiho et al. m- We fitted our model to data 
at 273.15K and 295.65K, to span the temperature range of CCS pipeline opera¬ 
tion. As the Fandiho et al. data do not have corresponding coexisting volumes, 
we required estimated volumes. Eor 273.15K we were able to validate our vol¬ 
ume estimates against the limited CO 2 -H 2 data from at this temperature, 
enabling us to successfully obtain H 2 parameters. However, the temperature 
change to 295.65K was sufficiently large that we were unable to successfully 
fit the Eandino et al. mol fraction data at this temperature using our volume 
estimate. Thus we allowed simulated annealing to adjust the parameters of our 
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e 

ao 

OCl 

OHa 

-0.919008 

0.910607 


-1.00994 

0.814616 

CHa 

-16.5057 

18.3776 

dHa 

0.0204381 

-0.0227559 

SHa 

-16.2801 

18.1264 

/Ha 

10.1075 

-10.2491 

5Ha 

-4.5342 

4.38617 


Table 5: Coefficients describing the temperature dependence of the H 2 model 
parameters in eqn (14). 


volume estimation technique in situe, during the fitting, to compensate for the 
lack of volume data. See |B.3| for details. We followed the same fitting procedure 
as the previous two cases and obtained the temperature dependence of each 
parameter via eqn (14). The coefficients of this linear temperature rule for each 
of the parameters are listed in table 

A comparison between the experiments and our model for the coexisting mol 
fraction at both temperatures is shown in figure]^ Also shown are the results 
from the GERG EoS. Our model captures these data very accurately, with 
the only disagreement being a slight discrepancy about the critical point for the 
higher temperature data. In contrast, the GERG EoS substantially over-predicts 
the critical point pressure at 295.65K, leading to large discrepancies that persist 
throughout the pressure range. Furthermore, at 273.15K the GERG EoS fails to 
capture the vapour mol fraction at moderate and higher pressures. Our model 
produces significantly closer agreement to these data than the GERG EoS. To 
our knowledge, the SAFT, PG-SAFT and EOS-GG have not been compared to 
these or comparable data. 


7 Discussion 

In this article we have developed a general framework for producing pressure- 
explicit EoS for impure CO 2 , aimed at modelling for GGS transport. Under our 
approach the mixture fugacity integral, required for coexistence calculations, 
can be computed from the pure fluid fugacity and mixing rules without fur¬ 
ther integration. This assists ongoing development of EoS in response to new 
data and computational requirements, as it allows convenient replacement and 
modification of terms in the EoS and mixing rules. 

We used our method to generalise a recent EoS for pure GO 2 to binary mix¬ 
tures. We used simulated annealing to fit model parameters to experimental 
data for mixtures with N 2 , O 2 and H 2 . We captured the coexistence data by 
imposing a term in the fitting that penalises differences in the predicted fugacity 
at the experimentally determined coexistence points. Where volumetric mea¬ 
surements were unavailable for coexistence, we developed, tested and exploited 
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Figure 6: Comparison of experiments and model predictions for the VLE 
coexisting mole fraction for CO 2 -H 2 mixtures. Experimental data from Fandino 
et al. [20] . Also shown are the results from the GERG EoS |0] ■ 


a volume estimation method. These estimated volumes allowed us to impose 
appropriate coexistence behaviour on our model. 

For CO 2 -N 2 mixtures, there was comprehensive and high quality literature 
data across virtually the entire CCS-pipeline window. Such data allowed us 
to achieve very good fitting with our model to both coexistence and pressure- 
volume data. For this system our model has comparable agreement to the 
GERG, but is slightly less accurate. However, for GO 2 -O 2 mixtures our model 
outperforms the GERG EoS, having more accurate coexisting mol fractions 
predictions at 273.15K and 288.15K. Our model is also more accurate than 
the SAFT EoS at 288.15K and has a similar level of agreement to the recent 
EOS-CG[T5|. We also successfully compared to coexisting mol fraction data 
for GO 2 -H 2 mixtures. Here our model is significantly more accurate than the 
GERG EoS, particularly at higher temperatures. 

A key issue for GGS modelling is the computational complexity of EoS. Some 
numerical codes use lookup tables of thermophysical properties, whereas other 
codes use direct, in-situe evaluation of EoS and so require very fast EoS. Our 
EoS has a moderate number of terms, with inexpensive linear relationships for 
the mixing rules and the temperature-dependence of the impurity parameters. 
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Furthermore, the pressure expression of our EoS contains only rational functions 
and so evaluation of this part of the EoS is computationally cheap. The required 
computational efficiency and whether or not to use lookup tables will depend 
upon the individual computational task. Ideally, a new EoS would be tailored to 
the computational and accuracy requirements of the application, a task towards 
which this work contributes. 

A persistent weakness for many EoS is modelling coexistence at pressures 
and temperatures in the vicinity of the critical point. Although our model is 
more accurate here than the GERG EoS for GO 2 -O 2 and GO 2 -H 2 , it has small 
discrepancies in this region for all impurities studied herein. This consistent 
behaviour suggests a systematic difficulty in our current approach, either in 
locating suitable parameters to describe the vicinity of the critical pressure or in 
the mathematical terms used in the EoS; and perhaps a difhculty in these factors 
across all current EoS. We discuss below, improved fitting algorithms that we are 
currently working on to address this problem. A further possibility to improve 
these issues is new mixing rules or modified terms in the EoS. We note that 
our EoS framework will make formulating these modifications straightforward 
and convenient. However, we leave exploration of these ideas to forthcoming 
studies. 

There are a number of weaknesses to our current approach. Firstly, there is 
a limit on the number of parameters that can be htted by simulated annealing to 
typical experimental data. To control the number of parameters in each search, 
we htted individual temperatures separately. Gonsequently, we could only ht 
temperatures where VLE mol fractions and volumes were available. Although 
this was mitigated somewhat by our volume estimation techniques, we still had 
to exclude some relevant data because the literature experiments were not sufh- 
ciently comprehensive at the relevant temperature. However, improved volume 
estimation or a modihed approach to htting coexistence data may improve the 
htting, especially for impurities where mixture data are sparse. 

Our approach opens up numerous possibilities to address the above weak¬ 
nesses. A key issue is the effectiveness of the parameter Htting algorithms. 
Effective non-linear optimisation is notoriously difficult as parameter-search al¬ 
gorithms often get trapped in deep local minima or struggle to explore effectively 
a complicated error function (such as eqn (HU))- In this work these factors lim¬ 
ited the number of parameters we could Ht simultaneously. These issues could 
be addressed by improved Htting algorithms. For example parallel tempering 
enables search algorithms to escape deep local minima [5^ and Riemann man¬ 
ifold, Langevin and Hamiltonian Monte Garlo methods are designed to cope 
with complicated error surfaces [55). Such improvements are likely to allow more 
model parameters to be Htted simultaneously. This will enable the temperature 
rules for model parameters to be Htted directly, thus meaning that experiments 
spread across many different temperature can be captured simultaneously. This 
may improve the EoS’s temperature-dependence around the critical point. 

There is also a central role for methods to address missing data or sparse 
volumetric data. We relied on our volume estimation techniques to stand-in for 
missing VLE data. A more accurate and physically-based estimate for these 
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volumes could be obtained from molecular simulation|561 15]. provided suitable 
molecular force fields describing the interaction of CO 2 with relevant impurities 
can be determined. It is also possible to treat missing VLE data as model 
parameters and learn these during the htting process, bypassing the need to 
explicitly estimate missing coexistence volumes. This introduces many new 
parameters to the fitting and so is contingent on solving the problem of fitting 
many parameters simultaneously, either by the techniques mentioned above, or 
otherwise. 

Reasonable solutions to the above issues will enable extensions of this work, 
targeting features that are useful to CCS applications. These include modelling 
a wider range of species of impurities, and ternary and higher order mixtures. 
Another future extension is uncertainty quantifications, which will address un¬ 
certainty due to sparsity of data for many impurities and guide the design of 
future experiments. 


8 Conclusions 

In this article we have developed a general framework for producing pressure- 
explicit EoS for impure CO 2 , aimed at modelling for CCS transport. This 
approach allows convenient modification of terms in the EoS and the mixing 
rules. We used our method to generalise a recent EoS for pure CO 2 to binary 
mixtures with N 2 , O 2 and H 2 , introducing 14 new model parameters per impu¬ 
rity. Our model pertains to pressures up to 16MPa and temperatures between 
273K and the critical temperature of pure CO 2 , Tc- For CO 2 -N 2 mixtures, our 
model has comparable agreement to the GERG, but is slightly less accurate. For 
CO 2 -O 2 mixtures our model outperforms the GERG EoS, for coexistence data 
at 273.15K and 288.15K. Here, our model is also more accurate than the SAFT 
EoS at 288.15K and has a similar level of agreement to the recent EOS-CG. We 
also successfully compared to coexisting mol fraction data for GO 2 -H 2 mixtures. 
Here our model is significantly more accurate than the GERG EoS, particularly 
at higher temperatures. We note that it is possible that some pipelines may op¬ 
erate above Tc, in some regions, which our pure GO 2 model does not currently 
account for. Fitting this region is more straightforward than the subcritical 
region due to the lack of coexistence. Furthermore, pipeline failure modelling 
requires an EoS that is accurate down to the GO 2 triple-point temperature, 
which is below the range we have explored here. Although our methodology is 
appropriate to these extensions, we leave them to future work. 
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A Mixture Fugacity in Equations of State 

We begin with the following fugacity expressions from Orby et aim-- equation 
(2.3.9) for the fugacity coefficient of a pure fluid, \n.(j), is 


In 


f*{T*,V* 


= \n(j) = 


RT* 


— -—dr*-\nZ + {Z-l), (15) 


where V*, N and Z = P*v*/RT* denote the volume, total number of moles and 
compressibility, respectively; and equation (2.3.1) for the fugacity coefficient of 
species i in a mixture is 


In 


ft 


p=¥ 


= In (/)i = 


1 




RT* 


RT* 

V'* 


dP* 


ON, 


dW*-lnZ, (16) 


* / T*,V,Nj^i 


where Ni and Xi are the total number of mols and mol fraction, respectively, of 
species i. All variables in dimensional units are notated with stars. 

Pressure-explicit EoS are of the form P* = P{v*,0*{Ni)), where v* = ^, 
so we can use the chain rule to write 


fdP* 


\dN, 




gp* y* 

'W*~Nfr 


^ dP* de* 


(17) 


where Np is the number of model parameters. Substituting into equation (161 
gives 


In^i = 


1 


rV* 


RT* 


RT* 


gp* y* \ 1 

dW*-lnZ - 


V* dV* Nt 


RT* 


^ de* 


rV* 


i=i 


dP* 

im 


dW 


( 18 ) 

Integration by parts shows that the term inside the brackets is equal to In (/> 
(by eqn ( |15[ )). Furthermore, defining F to be the integral from eqn (15), i.e. 
F = = luf/) — ((Z — 1) — In^), allows us to write. 


-N^ = ^ 
de* RT* 


dP* 

m* 


dV'- 


(19) 


Substituting (191 into equation (18) gives 




dF de* 

In (j)i = \n(j) + N ^ 


de* dN, 


( 20 ) 


The mixing rule will be of the form 0* = 0*{xi, {xky}), where Xi is the mol 
fraction of species i and the Xk are the mol fractions of all of the remaining 


20 






























species. Hence, we must now convert the derivative 8Ni in eqn ( |20[ ) into a 
derivative with respect to the mol fractions. We have Xi = — where 

Nother is the total number of mols of all species other than i. Similarly for 
k ^ i, we have Xk = -• Differentiating d* with respect to W gives 


Ni+N,th. 

dd* 

m ' 


do* dxi 
dxj dNj 


dOj dxk 


k^i 


dxk dNi' 


Differentiating the expressions for Xi and x^ allows us to write 

. dx 


N- 


dN,, 


= 1 — Xi and 


fi'T'h. 

dN,, 


Substituting equations (21) and (22) into equation (^20| gives 


, ? sr / a«; 




Ns: 

do* \ dxi ^ dxk 

0 = 1 3 \ ® k=l 


( 21 ) 


( 22 ) 


(23) 


We note here that, since the fugacity coefficients are dimensionless, then the 


mixing rules 6* are the only dimensional quantities in eqn (23). Furthermore 


because of the arrangement of the 6* terms in (23), the expression is unchanged 
by any constant scaling of the parameters 0*. Thus we conclude that the form 


of eqn (23) is unchanged by any non-dimensionalisation that applies a constant 


scaling to the model parameters. 


B Estimation Molar Volumes for Mixtures 

As noted previously, we observed that a lot of literature VLE data lacks measure¬ 
ment of the coexisting volumes. We thus developed from scratch a novel process 
to estimate both the liquid and vapour volumes in compensation for those data 
missing in the literature. This was primarily because our method of fitting 
the parameters required all elements of the thermodynamic description (tem¬ 
perature, pressure, volumes and compositions) to be present. We acknowledge 
that experimental measurements are clearly preferable to estimated volumes 
and have only used this technique when measurements are not available. 

B.l Estimating the Coexisting Liquid Volume 

We estimated the mixture liquid volume rtuq^Mix by taking a weighted average of 
the CO 2 volume at the pressure of interest and the pure impurity volume at the 
same pressure, weighted by the impurity concentration raised to an empirical 
power a. The pure data was available in all cases using standard reference 
EoS[8llS8l|59l|60j accessed from the NIST website [52]. Our empirical expression 
for the estimated volume of the coexisting mixture, uiiq^MiX: at the required 
temperature and pressure is, 

Wliq.MIX = (1 — 2;iiq)^C02 + a^UqUimp, (24) 
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where cciiq is the impurity concentration of the coexisting liquid, VCO 2 Wimp 
are the pure CO 2 and pure impurity molar volume at the temperature and pres¬ 
sure of interest, and a is an empirically determined constant, included to cali¬ 
brate this estimation for different mixtures, whose values for different mixtures 
are given in Table 

B.2 Estimating the Coexisting Vapour Volume 

The simple weighted average of eqn ( [24| ) is ineffective at estimating the coex¬ 
isting vapour volume, due to the higher impurity fraction and its influence on 
the CO 2 in the vapour phase. Thus a slightly more complicated expression is 
required. We used the following expression to estimate the coexisting vapour 
volume for the vapour,Uvap, mix 

V vnr^ -4- V ( '^^a.p.COs " t^liq.COa \ ^ ^ '|7^ 

^vap,MIX — ^vaD^C02 ' ^imp 1 ^vap/ ' V^vap ^liqj I ■> 

\ '^imp '^liq,C02 / 

(25) 

where ccvap is the impurity concentration of the coexisting vapour, Uvap,C 02 
and Wiiq,c 02 tde coexisting vapour and liquid volumes for pure CO 2 ,and 
/3 and 7 are two further empirically determined constants, given in Table 
The three empirical exponents depend only on the type of impurity and we 
determined these by fitting to the coexistence volume measurements used in 
this work. When fitting our EoS, we used experimental measurements for the 


Mixture 

a 

/3 

7 

C 02 “N 2 

1.25 

0.68 

1.01 

CO 2 -O 2 

1.55 

0.47 

0.73 

CO 2 —H 2 

1.51 

0.48 

1.19 


Table 6 : A summary of the values for the empirical exponents obtain by fitting 
to coexisting volume measurements.. 


coexisting volumes whenever available and used estimated volumes otherwise. 
The benefit of this method is that it allows a substitute value for the missing 
coexisting volumes, which often were not quoted, to be estimated. The technique 
requires only the pure density data, which was readily available through the 
NIST website [52] and the coexisting molar fractions, which are widely reported 
from experiments. A comparison between our estimation technique and co¬ 
existing volume measurements used in this work are shown in figures [^ and 

El 

B.3 Estimating the CO 2 -H 2 volumes close to the critical 
point 

For 273.15K we were able to validate our volume estimates against the limited 
CO 2 -H 2 data from ref [53| at this temperature, enabling us to successfully 
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obtain H 2 parameters (see figure]^. However, the temperature step to 295.65K 
was sufficiently large that we were unable to successfully fit the Fandino et al. 
VLE data at this temperature using our volume estimate. Thus we took a 
slightly more flexible form for the volume estimate 


i^vap 


.Mix(^’) = (fvap,C02 ~ ^'crit) 


Tcrit ~ P 

Pent — PcOq 




and 


(26) 


'Cliq,MIx(T’) = (^^liq.COa ~ ^crit) 


Tcrit ~ P 

J^crit - PC02 


"'^crit '^liq,C02 
^vap,C02 ^crit 


+ 'Ccrit, (27) 


where PCO 2 is the pure CO 2 coexistence pressure at 295.65K and Perit is the 
critical point pressure for CO 2 -H 2 mixtures at 295.65K, which Fandino et al. 
evaluated as PP^ = 14.655 MPa from their data. The parameters to be fitted 
are the exponent (j) and the mixture critical volume, ferit- We allowed simulated 
annealing to adjust the parameters of these two parameters, at the same time as 
the model parameters, during the fitting at 295.65K. This process gave values 
of ucrit = 0.320291 and (j) = 1.02222. 
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C02-N2 Data and Volume Estimation at 273.150K 



C02-N2 Data and Volume Estimation at 288.150K 
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Figure 7: Behaviour of the volume estimation method compared to literature 
data at 273.15K (top) [5S] and 288.15K (bottom) [32] for a CO 2 -N 2 mixture 
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C02—02 Data and Volume Estimation at 273.150K 
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Figure 8: Behaviour of the volume estimation method compared to literature 
data at 273.15K for a CO 2 -O 2 mixture 


C02-H2 Data and Volume Estimation at 273.150K 
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Figure 9: Behaviour of the volume estimation method compared to literature 
data at 273.15K for a CO 2 -H 2 mixture 
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